function y=Ryk(yk1, yk2, k, nt, N)
    y = zeros(2, 2);
    for m=0:nt-1
        y = y + [yk1(N*m + 1 + k);yk2(N*m + 1 + k)] * [yk1(N*m + 1 + k)' yk2(N*m + 1 + k)'];
    end
    y = y ./ nt;
end
